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Abstract:  This  paper  describes  the  optimization  of  a  load-bearing  thermal  in¬ 
sulation  system  characterized  by  hot  and  cold  surfaces  with  a  series  of  heat  in¬ 
tercepts  and  insulators  between  them.  The  optimization  problem  is  represented 
as  a  mixed  variable  programming  (MVP)  problem  with  nonlinear  constraints,  in 
which  the  objective  is  to  minimize  the  power  required  to  maintain  the  heat  in¬ 
tercepts  at  fixed  temperatures  so  that  one  surface  is  kept  sufficiently  cold.  MVP 
problems  are  more  general  than  mixed  integer  nonlinear  programming  (MINLP) 
problems  in  that  the  discrete  variables  are  categorical;  i.e.,  they  must  always  take 
on  values  from  a  predefined  enumerable  set  or  list.  Thus,  traditional  approaches 
that  use  branch  and  bound  techniques  cannot  be  applied. 

In  a  previous  paper,  a  linearly  constrained  version  of  this  problem  was  solved 
numerically  using  the  Audet-Dennis  generalized  pattern  search  (GPS)  method 
for  MVP  problems.  However,  this  algorithm  may  not  work  for  problems  with 
general  nonlinear  constraints.  A  new  algorithm  that  extends  that  of  Audet  and 
Dennis  by  incorporating  a  filter  to  handle  nonlinear  constraints  makes  it  possible 
to  solve  the  more  general  problem.  Additional  nonlinear  constraints  on  stress, 
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mass,  and  thermal  contraction  are  added  to  that  of  the  previous  work  in  an  effort 
to  find  a  more  realistic  feasible  design.  Several  computational  experiments  show  a 
substantial  improvement  in  power  required  to  maintain  the  system,  as  compared 
to  the  previous  literature.  The  addition  of  the  new  constraints  leads  to  a  very 
different  design  without  significantly  changing  the  power  required.  The  results 
demonstrate  that  the  new  algorithm  can  be  applied  to  a  very  broad  class  of  op¬ 
timization  problems,  for  which  no  previous  algorithm  with  provable  convergence 
results  could  be  applied. 


1  Introduction 

In  a  thermal  insulation  system,  heat  intercepts  are  often  used  to  minimize  the  heat  flow  from 
a  hot  to  a  cold  surface.  Figure  1  illustrates  an  example  of  such  a  system  of  fixed  length  L, 
in  which  power  is  applied  to  maintain  each  intercept  i  at  a  specified  cooling  temperature 
Ti,i  =  1,2, ...  ,n  (we  use  T  here  to  distinguish  it  from  the  general  temperature  variable  T 
used  later).  An  insulator  of  thickness  x,  is  placed  between  each  pair  of  intercepts  i  —  1  and  i, 
with  the  convention  that  i  —  0  and  i  =  n  + 1  represent  the  cold  and  hot  surfaces,  respectively, 
so  that  To  =  Tc  and  Tn+\  —  Th ■  Note  that  each  insulator  in  Figure  1  may  have  a  different 
cross-sectional  area.  The  design  variables  for  the  system  include  the  number  and  cooling 
temperatures  of  the  intercepts,  and  the  insulator  types  and  thicknesses.  Furthermore,  we 
assume  that  the  system  must  be  load-bearing,  meaning  that  the  insulators  act  as  mechanical 
supports;  thus,  only  solid  materials  can  be  used. 


Tn+ i  —  Th 

T  i+ 1 
Tt 

Ti-i 

T0  =  TC 


Figure  1:  Schematic  of  a  Thermal  Insulation  System 

Variations  of  this  type  of  system  used  in  cryogenic  engineering  applications,  such  as  su¬ 
perconducting  magnetic  energy  storage  systems  and  space  borne  magnets,  have  been  studied 
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by  several  authors.  Hilal  and  Boom  [22]  used  a  gradient-based  optimizer  to  minimize  power 
for  n  —  1,  2,  and  3  intercepts  with  two  choices  of  insulators  of  constant  cross-sectional  area, 
but  without  mixing  insulator  types  within  the  system.  Hilal  and  Eyssa  [23]  studied  the  same 
problem,  but  with  variable  cross  sectional  area  for  the  mechanical  supports.  In  considering 
systems  with  more  general  types  of  insulation,  Chato  and  Khodadidi  [14J  sought  to  minimize 
entropy,  similar  to  the  formulation  given  by  Bejan  [10] .  Other  related  attempts  to  optimize 
the  design  of  these  systems  are  found  in  [30J,  [35J,  and  [41J.  An  actual  implementation  of 
these  types  of  systems  for  the  Large  Hadron  Collider  (LHC)  project  is  discussed  in  three 
technical  reports  [18,  24,  31]. 

While  all  of  these  studies  vary  in  geometry  and  fidelity  of  the  underlying  models,  none 
of  them  optimizes  with  respect  to  the  number  of  intercepts  and  types  of  insulators.  These 
variables  are  referred  to  as  categorical,  meaning  that  they  must  take  on  values  only  from  a 
predefined  list  or  discrete  set  of  categories,  or  else  the  design  cannot  be  evaluated.  The  term 
appears  more  commonly  in  the  field  of  statistics,  in  which  each  observation  of  a  sample  is 
assigned  a  category  based  on  what  interval  its  value  lies  in  (e.g.,  see  [16]).  While  categorical 
variables  can  take  on  numerical  values,  the  values  may  not  have  any  real  meaning  or  inherent 
ordering  (e.g.,  1  =  steel,  2  =  aluminum,  etc.).  The  difficulty  these  variables  present  is  that 
traditional  methods  (most  notably,  any  branch  and  bound  approach)  for  solving  mixed  inte¬ 
ger  problems  cannot  be  applied  because  this  restriction  on  the  categorical  variables  cannot 
be  relaxed. 

Kokkolaras,  Audet,  and  Dennis  [26]  extended  the  Hilal  and  Boom  model  to  include  the 
number  of  intercepts  and  type  of  insulators  as  actual  variables  in  the  model,  thus  allowing 
a  mixture  of  different  insulator  types  in  the  system.  This  required  an  extension  of  the 
theory  [7]  to  be  able  to  deal  with  categorical  variables.  They  achieved  a  65%  reduction 
in  required  refrigeration  power  from  that  of  [22]  by  using  the  Audet- Dennis  mixed  variable 
generalized  pattern  search  (MGPS)  algorithm  [7]  to  solve  this  bound  constrained  mixed 
variable  programming  (MVP)  problem.  However,  other  than  choosing  the  specific  list  of 
possible  insulators,  they  did  not  consider  certain  load-bearing  aspects  of  the  problem,  such 
as  thermal  expansion,  system  mass,  and  stress,  because  these  are  modelled  as  nonlinear 
constraints. 

Since  the  MGPS  algorithm  does  not  handle  nonlinear  constraints,  a  new  GPS  algo¬ 
rithm  [2]  is  applied,  in  which  a  filter  is  added  to  the  MGPS  algorithm  to  handle  the  nonlinear 
constraints.  This  paper  describes  the  formulation  of  the  optimization  problem  (Section  2), 
the  extensions  that  have  been  made  to  the  Audet-Dennis  algorithm  (Section  3),  the  imple¬ 
mentation  of  the  algorithm  for  solving  the  problem  (Section  4),  and  computational  results 
(Section  5). 
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2  Basic  Model  of  Thermal  Insulation  Systems 

In  this  section,  we  reformulate  the  Kokkolaras  et  al.  model  (as  extended  from  that  of  [22J )  to 
include  constraints  on  thermal  expansion,  stress,  and  mass.  Much  of  this  discussion  comes 
directly  from  the  presentation  in  [26],  although  some  of  the  discussion  of  stress  is  similar  to 
that  of  [23] .  The  new  model  can  be  expressed  as 

min  f(n,I,x,T) 

(n,I,x,T)e  x  _  (1) 

subject  to  g(n,  I ,x,T)  <  0, 

with  the  following  nomenclature: 

•  n  is  the  number  of  heat  intercepts,  with  the  convention  that  the  cold  and  hot  walls  are 
numbered  0  and  n  +  1,  respectively. 

•  I  E  In+1  is  the  set  of  insulators  used,  where  represent  the  insulator  type  between 
intercepts  i  —  l  and  i,  and  1  denotes  the  finite  set  of  possible  insulator  types. 

•  x  G  3?n  is  the  vector  whose  2-th  component  is  the  thickness  of  the  i- th  insulator,  with 

the  convention  that  xn+i  =  L  —  i  xi- 

•  T  G  3?”  is  the  vector  whose  i-th  component  is  the  temperature  of  the  i-t\i  intercept, 
with  the  convention  that  T$  =  Tq  and  Tn+\  =Th- 

•  The  feasible  region  X  is  defined  by  the  following  linear  and  categorical  constraints: 


ti  G  {1,2,...,  7imax  { , 

(2) 

I  eln+\ 

(3) 

y ^Xi<  l, 

2=1 

(4) 

Xi  >  0,  i  —  1, . . . ,  n 

(5) 

Ti- 1  <Ti<  Ti+ 1,  2  =  1,.. 

• ,  n. 

(6) 

A  difficulty  in  solving  this  problem  is  that  the  dimension  of  the  vectors  /,  x,  and  T  depend 
on  the  variable  n.  For  any  value  of  n,  there  are  n  +  1  other  categorical  variables  and  2 n 
continuous  variables,  yielding  a  total  of  3n  +  2  variables. 

2.1  Objective  Function 

The  objective  function  represents  the  total  refrigeration  power  of  the  system;  thus, 

n 

/  =  EP- 

i=l 


(7) 
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where  Pi  is  the  power  applied  to  intercept  i,  i  =  1,2, ...  ,n.  The  required  power  to  keep 
intercept  i  at  a  fixed  temperature  Tt  is  given  by 

Pi  =  Ci  (qi+ 1  -  qi),  (8) 

where  Ci  (a  function  of  temperature  T)  is  the  thermodynamic  cycle  efficiency  coefficient  at 
intercept  i,  and  qi  represents  the  heat  flow  from  intercept  i  to  i  —  1. 

In  general,  heat  flow  q  through  a  volume  is  governed  by  Fourier’s  law, 


qdx  =  Ak(T)dT , 


(9) 


where  A  (a  function  of  the  spatial  coordinates  in  the  z-y  plane  perpendicular  to  the  x 
coordinate)  is  the  cross-sectional  area  of  the  volume,  and  k  (a  function  of  temperature)  is 
the  effective  thermal  conductivity  of  the  volume. 

For  the  problem  at  hand,  the  heat  flow  qi  from  intercept  i  to  i  —  1  is  given  by 

A-  fTi 

Qi  =  —  k(T;Ii)dT ,  i  =  1,  2, . . . , n  +  1,  (10) 

xi  Jp-i 

where  Ai  denotes  the  cross-sectional  area  of  insulator  i,  and  the  thermal  conductivity  k  is 
a  function  of  both  the  temperature  T  and  the  type  of  insulator  and  is  assumed  to  be 
isotropic. 

By  substituting  (10)  into  (8),  the  objective  function  can  be  expressed  by  (7),  with 


Pi 


Tlj+i 
%i-\- 1 


k(T-Ii+1)dT 


A  •  rT> 

—  /  k(T ;  Ii)dT 
xi  M. ! 


(11) 


2.2  Nonlinear  Constraints 

The  inequality  in  (1)  represents  the  new  nonlinear  constraints  added  to  the  MVP  problem 
of  |26j,  which  we  use  to  numerically  test  the  new  filter  MGPS  (FMGPS)  algorithm.  Although 
we  add  three  types  of  constraints;  namely  mass,  stress,  and  thermal  expansion,  we  model 
the  stress  constraint  as  an  implicit  equality  constraint  so  that  we  can  eliminate  insulator 
cross-sectional  areas  as  design  variables. 

The  first  constraint  concerns  the  overall  mass  of  the  system.  This  constraint  may  actually 
be  budgetary  in  nature,  as  larger  amounts  of  insulation  material  increase  the  overall  cost. 
Since  the  weight  of  each  insulator  can  be  expressed  as  a  product  of  the  density  of  the  insulator 
material  and  its  volume,  we  have  the  constraint, 

n 

J2Pi(Ii)AiXi<m 

max  i 

i= 1 


(12) 
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where  Pi(Ii)  represents  the  density  of  material  used  as  insulator  i,  and  mmax  is  the  maximum 
allowable  mass  of  the  system.  Note  that  setting  each  volume  to  AiXi  assumes  that  the  area 
of  the  insulator  is  constant  throughout  the  interval;  however,  this  is  not  an  unreasonable 
assumption. 

We  also  assume  that  the  system  must  be  capable  of  bearing  a  specified  load  F.  The  stress, 
defined  as  load  per  unit  area,  must  not  be  allowed  to  exceed  a  certain  level.  In  this  case,  we 
constrain  the  stress  applied  to  each  insulator  to  be  no  greater  than  the  tensile  yield  strength 
of  that  insulator  (thus,  we  assume  that  the  load  is  suspended  from  the  system,  rather  than 
resting  on  top  of  it).  The  tensile  yield  strength  at  insulator  i,  denoted  by  <7j(T;/j),  is  a 
function  of  both  the  insulator  type  and  temperature,  and  its  minimum  over  the  temperature 
range  of  the  insulator  cr*  must  act  as  the  constraint  limit  for  the  allowable  stress.  Therefore, 
we  have  constraints  of  the  form 

4-  <  cri  =  min  {cri(T;Ii)  :  Tj_i  <  T  <  T,},  i  =  1,  2, . . . ,  n  +  1.  (13) 


The  difficulty  with  the  constraints  given  in  (12)  and  (13)  is  that  they  treat  the  areas  Ai, 
i  =  1,  ...,n,  as  additional  design  variables.  However,  there  is  a  convenient  and  perfectly 
legitimate  way  around  this  problem.  First,  observe  that  the  direct  relationship  between  Ai 
and  Pi  in  (11)  means  that  decreasing  the  cross-sectional  area  of  any  insulator  reduces  the 
power  applied  to  the  corresponding  intercept,  which  is  exactly  the  goal  of  the  optimization 
problem.  Thus,  at  an  optimal  point,  each  Ai  should  be  as  small  as  possible.  Furthermore, 
as  each  Ai  is  made  smaller,  the  stress  constraint  in  (13)  becomes  binding  -  and  must  be 
so  at  optimality.  Therefore,  we  can  assume  that  (13)  holds  with  equality  and  make  the 
substitution  Ai  =  j-  to  eliminate  Ai  as  a  variable  in  (12)  and  in  the  objective  function.  This 
yields  a  stress-mass  constraint  of  the  form 


i=  1 


®max 

F 


(14) 


The  final  constraint  is  one  associated  with  the  system’s  thermal  expansion,  or  as  in  our 
cryogenic  application,  thermal  contraction.  In  addition  to  moving  the  heat  intercepts  out  of 
their  optimal  position,  thermal  contraction  causes  additional  stress  on  the  materials  and,  if 
excessive,  can  cause  other  difficulties,  such  as  deformations  in  the  material.  The  development 
of  this  constraint  as  presented  here  is  adapted  from  |9j. 

Since  different  insulators  at  different  temperatures  exhibit  different  contraction  behaviors, 
we  must  treat  thermal  contraction  of  each  insulator  separately  as  a  change  in  its  thickness. 
The  constraint  can  then  be  expressed  as  a  weighted  sum,  where  each  insulator’s  weight  is 
simply  its  thickness  divided  by  the  total  length  of  the  system;  i.e., 


(15) 
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where  5  is  a  limit  on  the  percent  total  contraction  of  the  system. 

For  insulator  i,  e(T;  /,;)  denotes  the  unit  thermal  contraction  (a  function  of  temperature) 
from  intercept  i  to  any  point  between  intercept  %  and  i  —  1  and  is  computed  by  the  formula, 

fT' 

e{T-It)  =  J i  XidT,  (16) 

where  T  G  [Ti_i,Ti],  and  A  j  is  the  linear  coefficient  of  thermal  expansion  for  insulator 
type  i.  Pre-computed  values  for  e(T;/j)  are  available  in  lookup  tables  for  a  wide  range 
of  temperatures  and  for  several  material  types  [9,  34J.  Total  thermal  contraction  for  an 
insulator  i  is  given  by 


Ax,  S7iyT-I,)k(T-l,)dT 


Xi  Sri,  MU 

thus,  the  nonlinear  thermal  expansion  constraint  is  given  by 


E 

2—1 


e(T]  Ii)k(T;  Ii)dT 


It!-!  k{T]  h)dT 


—  )  < 


5 

100' 


(17) 


(18) 


The  resulting  optimization  problem  can  now  be  expressed  by  the  objective  function  in 
(7)  and  (11),  modified  to  eliminate  the  Ai  variables,  together  with  the  linear  constraints 
defined  in  (2)-(6)  and  the  nonlinear  constraints  defined  by  (14)  and  (18);  namely, 


min_  P  =  F  y  Ci 

(n,I,x,T) 

v  '  n — \ 


subject  to 


1 


k(T-h+l)dT 


1 


(JiXi 


k(T ;  Ii)dT 


E 


StI,  e(T'>  Ii)k(T\  I,)dT 
Jl;  ,  k{T\I,)dT 


E««)-< 

“7  Or* 


®max 

F 


<  l, 

2=1 

Xi  >  0,  i  =  1, ...  ,7i 

Ti- 1  <  T,  <  Ti+ 1,  i  =  1, . . .  ,n, 

n  G  *(1,2,...,  nmax  }■ , 

I  G  ln+1. 
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3  Pattern  Search  Algorithms  for  Mixed  Variable  Op¬ 
timization 

Generalized  Pattern  Search  (GPS)  is  a  class  of  derivative-free  optimization  algorithms  that 
was  first  defined  and  analyzed  for  unconstrained  problems  by  Torczon  [40J,  and  extended 
to  problems  with  bound  [27J  and  linear  [28J  constraints  by  Lewis  and  Torczon.  Audet  and 
Dennis  |7j  extended  the  work  of  Lewis  and  Torczon  to  handle  MVP  problems  with  bound  con¬ 
straints.  This  algorithm  was  demonstrated  numerically  in  [26J  .  Audet  and  Dennis  provide 
additional  theoretical  results  under  a  hierarchy  of  more  general  smoothness  conditions  |8j. 
A  straightforward  extension  of  [7J  for  linearly  constrained  MVP  problems  under  the  relaxed 
smoothness  assumptions  is  described  in  [2j. 

For  problems  with  nonlinear  constraints,  Lewis  and  Torczon  [29J  describe  an  augmented 
Lagrangian  GPS  approach  in  which  GPS  is  used  to  generate  approximate  solutions  to  a 
sequence  of  bound  constrained  augmented  Lagrangian  subproblems  (see  [15J ) .  As  an  al¬ 
ternative,  Audet  and  Dennis  |6J  introduce  a  filter  GPS  algorithm.  The  latter  cannot  quite 
guarantee  convergence  to  a  first-order  stationary  point  as  the  former  can,  but  it  applies  di¬ 
rectly,  rather  than  sequentially,  to  more  general  problems,  and  it  avoids  the  use  of  penalty 
parameters  and  Lagrange  multipliers. 

The  new  FMGPS  algorithm  used  here  is  an  extension  of  the  mixed  variable  GPS  algo¬ 
rithm  of  Audet  and  Dennis,  in  which  a  filter  is  added  to  handle  the  nonlinear  constraints. 
It  represents  a  generalization  of  all  the  work  cited  here,  and  is  described  in  detail  in  [4J  and 
[2]- 

3.1  Filters 

Fletcher  and  Leyffer  [19J  developed  the  filter  algorithm  as  a  way  to  globalize  sequential 
quadratic  programming  (SQP)  and  sequential  linear  programming  (SLP)  methods  without 
the  need  for  a  merit  or  penalty  function,  which  would  require  the  user  to  specify  the  relative 
weighting  of  optimality  versus  feasibility.  Proofs  of  convergence  to  a  first-order  stationary 
point  are  given  in  [20J  for  SQP  and  [21j  for  SLP  methods. 

In  filter  algorithms,  the  goal  is  to  minimize  two  functions,  the  objective  /  and  a  continuous 
aggregate  constraint  violation  function  h  that  satisfies  h(x)  >  0  with  h(x)  =  0  if  and  only 
if  x  is  feasible.  The  function  h  is  often  set  to  h(x)  =  ||G(a;)+||,  where  ||  -  ||  is  a  vector  norm 
and  C(x)+  is  the  vector  of  constraint  violations  at  x ;  i.e.,  for  i  —  1,  2, . . . ,  m,  Ci(x)+  =  Cfix) 
if  Ci(x)  >  0;  otherwise,  Ci(x)+  =  0.  In  particular,  if  the  squared  2-norm  is  used,  then  h 
inherits  whatever  smoothness  properties  C  possesses  |6J. 

Simply  defined,  a  filter,  denoted  T ,  is  a  set  of  trial  points  such  that  no  point  dominates 
any  other  in  the  set  with  respect  to  its  objective  and  constraint  violation  function  values.  In 
other  words,  given  any  two  points  x  and  y  in  the  filter,  either  f(x)  <  f(y)  or  h(x)  <  h(y), 
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but  not  both. 

In  constructing  a  filter  for  GPS,  we  put  two  additional  restrictions  on  T .  First,  we  set  a 
bound  on  aggregate  constraint  violation,  hmax,  so  that  each  point  y  <G  JF  satisfies  h(y)  <  /rmax. 
Second,  we  include  only  infeasible  points  in  the  filter  and  track  feasible  points  separately. 
This  is  done  in  order  to  avoid  a  problem  with  what  Fletcher  and  Leyffer  [19J  refer  to  as 
“blocking  entries”,  in  which  a  feasible  filter  point  with  lower  function  value  than  a  nearby 
local  minimum  prevents  convergence  to  both  that  minimum  and  a  global  minimum.  Tracking 
feasible  points  outside  of  the  filter  circumvents  this  uncommon  but  plausible  scenario.  With 
these  modifications,  we  now  refer  to  a  point  y  as  filtered  if  it  is  dominated  by  any  point  in 
the  filter,  satisfies  h(y)  >  hmax,  or  is  feasible  and  has  objective  function  value  greater  than 
the  incumbent  best  feasible  point.  A  point  that  is  not  filtered  is  referred  to  as  unfiltered. 

3.2  The  Filter  MGPS  Algorithm 

A  basic  pattern  search  algorithm  is  characterized  by  two  phases  per  iteration  -  a  global 
SEARCH  and  a  local  POLL,  in  which  trial  points  lying  on  a  carefully  constructed  mesh  are 
evaluated.  The  goal  of  the  search  is  to  adequately  sample  the  variable  space,  so  as  to 
quickly  identify  a  promising  region  containing  a  good  local  minimizer.  In  this  step,  any 
finite  strategy  may  be  employed  (including  none)  to  identify  trial  points,  as  long  as  they 
lie  on  the  mesh.  The  gives  the  user  great  flexibility  in  choosing  how  to  select  points.  For 
example,  any  of  the  following  strategies  may  be  used  in  the  search  step  to  identify  trial 
points  to  evaluate. 

•  Randomly  select  mesh  points  using  a  Latin  hypercube  search  [33,  38,  39J  or  orthogonal 
arrays  [36J; 

•  Use  a  popular  heuristic,  such  as  a  few  iterations  (be.,  generations)  of  a  genetic  algo¬ 
rithm; 

•  Optimize  a  less  expensive  surrogate  function  and  map  the  resulting  numerical  solution 
to  its  closest  mesh  point  (see  [12,  13J  as  examples). 

The  latter  option  is  popular  in  practice  for  optimization  problems  whose  objective  and 
constraint  functions  are  expensive  to  evaluate.  An  initial  investment  of  function  evaluations 
is  needed  to  form  good  surrogate  objective  and  constraint  functions,  but  then  optimization 
of  the  surrogate  problem  may  be  quickly  achieved  with  little  expense,  and  the  resulting  point 
mapped  to  the  mesh  should  yield  an  improved  design.  If  an  improved  design  is  not  found, 
the  surrogate  functions  are  updated  and  then  re-optimized  during  the  next  iteration.  Several 
studies  have  shown  success  using  this  approach  [5,  11,  12,  13J. 

If  the  search  step  does  not  produce  an  improved  mesh  point  (be.,  a  point  with  a  lower 
objective  function  value  than  the  current  iterate),  the  poll  step  is  executed,  in  which  the 
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mesh  neighbors  of  the  current  iterate  (called  the  poll  set)  are  evaluated.  This  step  is  per¬ 
formed  so  that  convergence  of  a  subsequence  of  iterates  is  guaranteed  to  produce  a  limit 
point  satisfying  certain  necessary  conditions  for  optimality.  The  careful  construction  of  the 
mesh  enables  us  to  converge  to  such  a  point.  The  poll  step  can  be  terminated  early  if  an 
improved  mesh  point  is  found. 

For  filter  GPS  algorithms,  the  mesh  and  poll  set  are  defined  in  terms  of  a  poll  center,  as 
opposed  to  the  current  iterate.  That  is,  at  each  iteration  k,  the  poll  center  pk  G  {pk,p{}  is 
chosen  as  either  the  incumbent  best  feasible  point  pk  or  the  incumbent  least  infeasible  point 
pTk  (/  and  F  denote  infeasible  and  feasible,  respectively).  Instead  of  seeking  a  trial  point 
with  a  better  objective  function  value  than  the  incumbent,  the  filter  GPS  algorithm  searches 
for  an  nn filtered  point  that  can  be  added  to  the  filter  -  ideally,  an  improved  incumbent  best 
feasible  or  least  infeasible  point. 

For  mixed  variable  problems,  each  iterate  Xk  =  (xk,  xk),  where  xrk  is  the  vector  of  continu¬ 
ous  variables  and  xk  is  the  vector  of  discrete  or  categorical  variables.  The  mesh  is  constructed 
as  the  direct  product  of  the  discrete  variable  space  Xd  with  the  union  of  a  finite  number  of 
lattices  in  the  continuous  variable  space  Xc,  but  translated  from  the  poll  center;  i.e., 

^max 

Mk  =  Xd  x  |J  {p%  +  A kDlz  G  Xc  :  z  G  (19) 

i=l 

where  the  mesh  size  parameter  Ak  is  a  positive  real  number  that  controls  the  fineness  of  the 
mesh,  and  Dl  denotes  an  n  x  nDi  (where  >  n )  matrix  whose  columns  positively  span 
the  continuous  variable  space.  Typically,  these  are  chosen  as  [/,—/]  or  [/,  — e],  where  /  is 
the  identity  matrix  and  e  is  the  vector  of  ones. 

The  poll  set  is  the  union  of  the  continuous  mesh  neighbors  P(pk )  of  pk  with  the  user- 
defined  set  of  discrete  neighbor  points  A f(pk)  of  pk.  The  continuous  mesh  neighbors  can  be 
expressed  as 


PM  =  {Pi  +  A kd  G  Xc  :  d  G  Dl},  (20) 

where  Dlk  C  Dl  is  the  set  of  poll  directions  at  iteration  k.  For  convenience,  we  use  the 
notation,  hTk  =  h^p}.)  >  0,  fk  —  f(p{),  and  fk  =  f{pk).  The  flexibility  in  selecting  either  of 
the  two  poll  centers  does  not  affect  the  convergence  theory  [3,  6j,  but  it  can  cause  convergence 
to  a  different  limit  point.  Audet  and  Dennis  |6j  (see  Example  7.1  there)  show  that  choosing 
one  or  the  other  exclusively  in  can  result  in  unfavorable  results  in  practice.  They  suggest 
that  one  good  strategy  may  be  to  alternate  between  the  two  poll  centers  each  time  a  mesh 
isolated  filter  point  is  found. 

We  should  note  that  polling  around  other  points  in  the  filter  is  certainly  allowed;  however, 
other  points  do  not  possess  the  same  convergence  properties  as  pk  and  pk  in  the  limit. 
Therefore,  any  such  polling  should  be  regarded  as  part  of  the  SEARCH  step,  so  that  the 
convergence  theory  is  preserved. 
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In  the  MGPS  method  of  Audet  and  Dennis  [7J,  if  the  poll  step  fails  to  yield  an  improved 
mesh  point,  an  EXTENDED  POLL  step  is  invoked  around  any  discrete  neighbor  whose  objective 
function  value  is  sufficiently  close  to  that  of  the  current  iterate  (he.,  “almost”  an  improved 
mesh  point).  This  is  done  with  the  goal  of  generating  an  improved  mesh  point. 

With  the  addition  of  nonlinear  constraints  to  the  problem,  the  filter  approach  seeks  a 
point  that  improves  either  the  objective  function  /  or  the  constraint  violation  function  h. 
Given  the  current  poll  center  p};  and  user-specified  extended  poll  triggers  >  £  >  0  and 
>  £  >  0  for  /  and  h,  respectively  (for  some  positive  constant  £),  we  perform  an  EXTENDED 
POLL  step  around  any  feasible  discrete  neighbor  whose  objective  function  value  is  within 
of  that  of  the  best  feasible  point,  or  around  any  infeasible  discrete  neighbor  whose  constraint 
violation  function  value  is  within  ££  of  that  of  the  least  infeasible  point  (without  exceeding 

An  ax  )  • 

The  bi-loss  graph  in  Figure  2  depicts  a  filter  with  its  possible  poll  centers  for  the  next 
iteration.  The  best  feasible  and  least  infeasible  points  are  indicated,  and  the  feasible  solutions 
lie  on  the  vertical  axis  (labelled  /).  The  dashed  lines  have  been  added  to  indicate  the  areas  for 
which  an  extended  poll  step  is  triggered.  If  a  feasible  discrete  neighbor  has  an  objective 
function  value  higher  on  the  axis  than  the  current  incumbent,  but  lower  than  the  horizontal 
dashed  line,  an  EXTENDED  POLL  step  is  performed  around  this  discrete  neighbor.  Similarly, 
an  EXTENDED  POLL  step  is  performed  if  a  filtered  discrete  neighbor  lies  to  the  right  of  the 
current  least  infeasible  solution,  but  left  of  the  vertical  dashed  line. 
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Mesh  isolated  filter  points:  Tk  C  Tk 


Figure  2:  A  Filter  for  GPS  with  Extended  Poll  Triggers. 

Similar  to  the  MGPS  algorithm  of  [7J,  the  extended  poll  step  generates  a  finite 
sequence  of  extended  poll  centers  during  a  single  iteration,  beginning  with  the  discrete 
neighbor  z/k-  However,  the  progress  of  this  step  is  also  different  from  that  of  [7j,  since  we 
seek  an  unfiltered  point,  as  opposed  to  a  simple  decrease  in  the  objective  function. 

At  each  step  in  the  EXTENDED  POLL,  there  are  three  possible  outcomes.  First,  if  the 
EXTENDED  POLL  yields  an  un filtered  point,  then  the  point  is  treated  exactly  the  same  as 
any  other  unfiltered  point;  namely,  the  point  is  added  to  the  filter,  the  filter  is  updated,  and 
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the  mesh  is  coarsened  or  left  unchanged.  Second,  if  the  extended  poll  fails  to  find  a  point 
with  a  better  objective  or  constraint  violation  function  value  than  the  current  extended  poll 
center,  then  the  EXTENDED  POLL  step  around  the  current  discrete  neighbor  is  terminated, 
and  extended  polling  begins  around  the  next  qualified  discrete  neighbor.  Once  all  such 
points  have  been  tested  with  the  same  result,  the  EXTENDED  POLL  step  is  terminated,  and 
the  mesh  is  refined. 

However,  the  third  possible  outcome  requires  a  new  construction;  that  is,  if  extended 
polling  around  a  discrete  neighbor  fails  to  yield  a  unfiltered  point,  but  improves  the  objec¬ 
tive  or  constraint  violation  function  value,  relative  to  the  discrete  neighbor  point ,  then  the 
extended  poll  step  should  be  allowed  to  continue,  since  finding  an  unhltered  point  is  still 
promising.  This  is  done  by  establishing  a  temporary  local  filter.  At  iteration  k ,  for  each 
discrete  neighbor  yk,  a  local  filter  iFfi(yk)  is  simply  a  filter  relative  to  the  current  EXTENDED 
POLL  step.  It  is  populated  initially  with  only  the  point  tjk  and  with  h^ax  =  min(h[.+^,  hmax). 
The  extended  poll  step  then  generates  a  finite  sequence  of  extended  poll  centers,  where 
each  is  chosen  either  as  the  best  feasible  or  least  infeasible  point,  relative  to  the  local  filter. 
Extended  polling  with  respect  to  yk  proceeds,  with  points  being  added  to  the  local  filter  as 
appropriate,  until  no  more  unhltered  mesh  points  can  be  found  with  respect  to  the  new  local 
filter,  or  until  an  unhltered  point  is  found  with  respect  to  the  main  hlter.  Once  either  of  these 
conditions  is  satished  (which  is  guaranteed  to  occur  for  a  hxed  mesh  size),  the  extended 
poll  step  ends,  and  the  main  hlter  is  appropriately  updated  with  the  points  of  the  local 
hlter,  which  is  then  discarded.  The  mesh  size  parameter  A k,  which  has  been  kept  constant 
throughout  the  step,  is  then  updated,  depending  on  the  success  of  the  SEARCH,  POLL,  and 
EXTENDED  POLL  steps  in  finding  an  unhltered  point  with  respect  to  the  main  hlter. 

If  a  trial  point  [i.e.,  any  point  evaluated  in  the  three  phases)  is  dominated  by  any  point 
in  the  hlter,  it  is  said  to  be  “filtered” .  If  all  the  trial  points  are  hltered,  then  the  current  poll 
centers  are  retained  and  the  mesh  is  refined  by  setting  A/j+1  <  A*,,  typically  A *.+1  =  A^/2. 
Otherwise,  the  trial  point  is  added  to  the  hlter,  the  hlter  is  updated  to  remove  any  points  that 
are  now  dominated  by  the  new  point,  and  the  mesh  is  coarsened  by  setting  A^+i  >  A k-  Note 
that  this  includes  the  choice  to  leave  the  mesh  unchanged.  More  precise  rules  that  govern 
how  A k  must  be  refined  or  coarsened  to  maintain  convergence  properties  of  the  algorithm 
are  given  in  |40j,  [8j,  and  [7j. 

The  EXTENDED  POLL  step  is  summarized  by  the  algorithm  in  Figure  3,  and  the  FMGPS 
Algorithm  is  summarized  in  Figure  4 


4  Computational  Model 


This  section  further  describes  the  optimization  problem  in  terms  of  modelling  decisions, 
material  data  sources,  and  other  problem  setup  issues. 


May  5,  2003 


13 


Extended  Poll  Step  at  Iteration  k 

Input:  Current  poll  center  pk,  filter  Pk,  and  extended  poll  triggers  and 

For  each  discrete  neighbor  yk  satisfying  the  EXTENDED  POLL  criteria,  do  the  follow¬ 
ing: 

•  Initialize  local  Elter  Tk  with  yk  and  h^ax  =  min  ihk  +  Set  yl  =  yk. 

•  For  j  =  0,1,2,... 

1.  Evaluate  /  and  h  at  points  in  Pk(y3k)  until  a  point  w  is  found  that  is 
unfiltered  with  respect  to  Fk  ,  or  until  done. 

2.  If  no  point  w  G  Pk(y{)  is  unfiltered  with  respect  to  JF^,  then  go  to  Next. 

3.  If  a  point  w  is  unfiltered  with  respect  to  Pk,  set  xk+i  =  w  and  Quit. 

4.  If  w  is  filtered  with  respect  to  Pk,  but  unfiltered  with  respect  to  JF^,  then 
update  Pk  to  include  w ,  and  compute  new  extended  poll  center  y{+1. 

•  Next:  Discard  JF^  and  process  next  yk. 

Figure  3:  Extended  Poll  Step  for  the  FMGPS  Algorithm 

4.1  Material  Data 

The  types  of  insulators  were  chosen  as  the  same  as  in  [26J;  namely,  nylon,  teflon,  fiberglass 
epoxy  (both  normal  and  plane),  6063-T5  aluminum,  1020  low-carbon  steel,  and  304  stainless 
steel.  For  each  of  these  materials,  a  substantial  amount  of  engineering  data  was  required. 
Thermal  conductivity  and  contraction  data  were  obtained  from  lookup  tables  in  |9J  and  [34J , 
while  material  densities  were  found  in  |37J  and  |34j,  and  tensile  yield  strength  data  were 
obtained  from  [17J  and  |32j. 

Thermodynamic  cycle  efficiency  coefficients  C$,  i  —  1,  2, . . . ,  n  (see  (8))  are  dependent  on 
temperature  as  follows: 

(  5,  if  T  <  4.2 K 

Ci  =  <  4,  if  4.2 K  <  T  <  71 K  (21) 

{  2.5,  if  T  >  71 K. 

In  order  to  make  the  model  as  accurate  and  efficient  as  possible,  cubic  splines  were 
used  to  fit  all  of  the  data  found  in  lookup  tables,  including  thermal  conductivity,  thermal 
contraction,  and  tensile  strength  data.  Numerical  integrations  were  performed  by  applying 
a  composite  Simpson’s  Rule,  with  nodes  matching  those  of  the  cubic  spline.  This  eliminates 
truncation  error,  since  Simpson’s  Rule  is  exact  for  cubic  polynomials  [25J. 
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Filter  Mixed  Variable  Generalized  Pattern  Search  -  FMGPS 

Initialization:  Let  Xo  be  an  undominated  point  of  a  set  of  initial  solutions.  Include 
all  these  points  in  the  filter  JF0,  with  hmax  >  h(x 0).  Fix  £  >  0  and  A0  >  0. 

For  k  =  0,1,2, ,  perform  the  following: 

1.  Update  poll  center  pk  G  {Pk-,p{}i  extended  poll  triggers  >  £  and  ££  >  £. 

2.  Compute  incumbent  values  =  /(pf ),  h[  =  h(p{),  f[  =  f(p{). 

3.  Search  step:  Employ  some  finite  strategy  seeking  an  unfiltered  mesh  point 

4.  POLL  step:  If  the  SEARCH  step  did  not  End  an  un filtered  point,  evaluate  / 
and  h  at  points  in  the  poll  set  Pk(Pk )  U  Af(pk)  until  an  unfiltered  mesh  point 
Xk+i  is  found,  or  until  done. 

5.  Extended  Poll  step:  If  search  and  poll  did  not  find  an  unfiltered  point, 
execute  the  algorithm  in  Figure  3  to  continue  looking  for  an  unfiltered  point 

%k+l  • 

6.  Update:  If  search,  poll,  or  extended  poll  finds  an  unfiltered  point, 

Update  filter  Tk+i  with  Xk+ 1,  and  set  A^+i  >  A&; 

Otherwise,  set  J-k+i  =  Pk,  and  set  A^+i  <  A*,. 

Figure  4:  FMGPS  Algorithm 

4.2  Choosing  Discrete  Neighbors 

The  neighborhood  structure  that  the  user  chooses  to  incorporate  essentially  determines  the 
definition  of  a  minimizer.  That  is,  when  the  solution  of  an  MVP  problem  is  found,  it  is  with 
respect  to  the  user-specified  discrete  set  of  neighbors.  If  a  neighborhood  structure  is  chosen 
so  that  all  other  sets  of  discrete  variable  values  are  neighbors,  a  more  global  solution  can 
be  obtained,  but  often  at  extraordinary  computational  cost.  On  the  other  hand,  severely 
restricting  the  size  of  the  set  of  neighbors  will  save  significant  computational  cost,  but  a  local 
optimizer  with  a  higher  objective  function  value  is  likely  to  be  obtained. 

In  order  to  make  proper  comparisons,  the  set  of  neighbors  we  use  for  this  problem  is 
exactly  the  same  as  was  used  by  Kokkolaras  et  al.  |26j.  It  includes  designs  in  which  the 
following  occur: 

•  The  type  of  insulator  between  any  two  heat  intercepts  is  changed  to  any  other  type, 
while  insulator  thicknesses  and  intercept  temperatures  remain  constant. 
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•  An  intercept  and  the  insulator  above  it  are  removed,  while  thicknesses  of  the  remaining 
insulators  are  increased  proportionally  (rounded  to  the  nearest  integer  multiple  of  the 
current  mesh  size)  to  fill  the  remaining  space. 

•  A  new  intercept  and  an  insulator  underneath  it  are  added  with  the  following  properties: 

—  The  type  of  insulator  is  the  same  as  the  one  below  it, 

—  The  cooling  temperature  is  set  to  the  average  of  the  two  intercepts  adjacent  to  it, 
rounded  to  the  nearest  integer  multiple  of  the  current  mesh  size, 

—  The  thickness  of  both  the  new  insulator  and  the  insulator  below  it  are  both  set  to 
half  of  that  of  latter,  rounded  to  the  nearest  integer  multiple  of  the  current  mesh 
size. 

Note  that  rounding  to  the  nearest  integer  multiple  of  the  current  mesh  size  is  necessary 
to  ensure  that  the  trial  point  lies  on  the  mesh. 


5  Computational  Results 

The  FMGPS  algorithm  described  in  Section  3.2  has  been  implemented  in  a  Matlab®  code, 
called  NOMADm  [1],  and  applied  to  this  problem.  We  now  present  results  of  several  NO- 
MADm  runs  and  show  that  the  approach  significantly  improves  the  design  of  Hilal  and 
Eyssa  |23j ,  and  is  comparable  to  the  results  of  Kokkolaras  et  al.  |26j ,  even  though  the  prob¬ 
lem  takes  on  the  additional  nonlinear  load-bearing  constraints. 

Table  1  shows  the  data  parameters  chosen  for  the  results  that  follow.  The  first  four  have 
values  identical  to  the  choices  of  Kokkolaras  et  al.  [26J ,  while  the  remaining  parameters  are 
unique  to  this  problem. 


Table  1:  Thermal  Insulation  Problem:  Model  Parameters 


Parameter 

Symbol 

Value 

Hot  surface  temperature 

Th 

300  K 

Cold  surface  temperature 

Tc 

4.2  K 

System  total  length 

L 

100  cm 

Maximum  number  of  intercepts 

^max 

10 

Load  placed  on  the  system 

F 

250  kN 

Maximum  total  system  mass 

®max 

10  kg 

Maximum  system  thermal  contraction 

5 

5% 
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To  match  the  setup  of  [26]  as  much  as  possible,  runs  were  performed  with  an  initial  mesh 
size  of  A0  =  10  and  terminated  when  the  condition  A*,  <  .15625  was  achieved.  The  mesh 
refinement  strategy  used  by  |26j  could  not  be  duplicated  by  NOMADm;  thus,  in  our  runs, 
we  refine  the  mesh  by  simply  dividing  the  mesh  size  parameter  Afc  in  half.  Coarsening  of 
the  mesh  was  not  performed. 

Extended  poll  triggers  for  the  objective  and  constraint  violation  function  were  set  at  one 
and  five  percent,  respectively,  the  former  being  consistent  with  [26J.  When  the  filter  logic 
of  the  FMGPS  is  applied  for  the  nonlinear  constraints,  polling  is  performed  around  the  best 
feasible  point. 

Also  consistent  with  [26],  no  SEARCH  is  used,  and  the  initial  design  consisted  of  one 
intercept  placed  exactly  in  the  middle  of  the  system  and  set  at  150  K,  with  a  nylon  insulator 
on  the  cold  side  and  a  teflon  insulator  on  the  hot  side. 

5.1  Validation 

The  software  and  function  hies  were  validated  by  mimicking  the  designs  of  [22]  and  [23] , 
running  these  problems,  and  comparing  the  designs.  In  both  of  these  previous  papers,  the 
authors  applied  their  optimizer  to  cases  of  1-3  heat  intercepts  with  insulators  made  of  either 
304  stainless  steel  or  plane-cloth  fiberglass  epoxy.  For  stainless  steel,  our  results  matched 
theirs  almost  exactly.  For  fiberglass  epoxy,  there  were  some  slight  differences  when  more 
than  one  intercept  is  used,  but  these  were  noted  in  [26J  as  well.  These  differences  were  most 
likely  caused  by  different  methods  in  computing  the  objective  function  integrals.  As  in  [26], 
we  used  cubic  splines  to  fit  thermal  conductivity  data  that  was  available  only  in  tabular 
form  for  specific  temperatures.  However,  rather  than  apply  a  Matlab  integration  routine, 
we  applied  our  own  implementation  of  Simpson’s  rule,  which  is  exact  for  cubic  polynomials. 
The  inaccuracy  is  more  visible  in  the  epoxy  results  because  thermal  conductivity  data  was 
only  available  at  four  temperatures,  as  opposed  to  the  18  different  temperatures  available 
for  stainless  steel. 

When  we  recomputed  the  results  of  [26]  to  validate  our  mixed  variable  logic,  we  converged 
to  a  different  design  with  a  similar  low  objective  function  value.  Table  2  shows  the  differences 
between  the  two  runs,  where  the  materials  cited  there  are  abbreviated  by  the  following:  N  = 
nylon,  E  =  epoxy  (normal),  Ep  =  epoxy  (plane),  and  T  =  teflon.  Note  that,  although  power 
is  optimized  in  our  software,  we  report  normalized  power  at  termination,  in  which  power  is 
multiplied  by  the  system  length  L  and  divided  by  the  smallest  cross-sectional  area  of  any 
insulator.  Previous  authors  have  expressed  results  this  way  so  that  designs  can  be  compared, 
independent  of  these  two  parameters.  We  keep  this  convention  for  the  same  purpose. 

While  the  numerical  integration  issue  just  described  can  lead  to  small  deviations,  we 
suspect  that  the  difference  in  mesh  refinement  strategies  led  to  a  different  local  optimizer 
along  a  different  path.  In  [26],  with  a  starting  mesh  size  of  Ao  =  10,  the  mesh  refinement 
strategy  was  to  divide  the  current  mesh  size  by  2l,  where  £  =  1,  2, ...  is  incremented  each 
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time  the  mesh  is  refined.  Since  NOMADm  is  currently  incapable  of  incrementing  £,  our  mesh 
refinement  consists  of  simply  dividing  the  current  mesh  size  by  two  (he.,  I  =  1). 


Table  2:  Thermal  Insulation  Problem:  MVP  Validation 


Problem: 

Original 

FMGPS  ReRun 

Power  (^): 

25.294  W/crn 

25.589  W/cm 

Insulators  Used: 

NNNNNNNEEET 

NNNNNNTEETT 

i 

Xi  (cm) 

Ti  (K) 

Xi  (cm) 

Ti  (K) 

1 

0.3125 

4.2188 

4.5313 

6.125 

2 

5.4688 

7.3438 

6.7188 

10.55 

3 

3.9062 

10 

4.8437 

14.35 

4 

6.5625 

15 

4.2188 

17.994 

5 

5.7812 

20 

7.3438 

24.969 

6 

5.1562 

25 

9.8438 

36.006 

7 

13.2812 

40 

24.948 

71.094 

8 

21.4062 

71.0938 

12.135 

116.88 

9 

8.5938 

101.25 

7.5 

156.88 

10 

9.2188 

146.25 

6.4063 

198.44 

11 

20.3125 

11.5105 

5.2  Adding  the  Nonlinear  Constraints 

The  nonlinear  constraints  were  added  to  the  runs  in  three  steps.  First,  we  simply  tested  the 
two  designs  from  Table  2  versus  the  new  nonlinear  constraints  and  found  that  the  thermal 
contraction  constraint  for  both  designs  was  violated  by  approximately  8%  (he.,  the  observed 
thermal  contraction  of  5.4%  exceeded  the  5%  threshold  by  8%).  This  suggests  that  a  new 
design  having  a  different  material  configuration  should  be  expected  as  the  new  constraints 
are  incorporated.  Second,  the  implicit  constraint  on  stress  (given  by  (13)  with  equality),  was 
added  to  allow  variable  cross-sectional  areas,  and  thus  match  the  formulation  of  Hilal  and 
Eyssa  [23J.  We  refer  to  this  as  the  partial  model.  Finally,  the  mass  and  thermal  contraction 
constraints  were  added  to  complete  the  full  model.  By  doing  so,  the  resulting  change  in 
required  power  represents  the  cost  of  satisfying  the  additional  load-bearing  constraints. 

Table  3  shows  the  results  for  the  partial  and  full  models  (columns  3  and  4),  along  with 
the  design  found  by  Hilal  and  Eyssa  |23j  (column  2).  For  each  run,  the  minimal  required 
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power  is  normalized  by  multiplying  by  the  total  system  length  and  dividing  by  the  smallest 
cross-sectional  area  of  an  insulator.  Following  the  power  and  the  material  configuration  for 
each  design,  insulator  thicknesses  and  heat  intercept  temperature  settings  are  listed. 


Table  3:  Thermal  Insulation  Problem:  Results 


Problem: 

Hilal  &  Eyssa 

Partial  Model 

Full  Model 

Power  (^A): 

53.2  W/cm 

24.551  W/cm 

23.768  W/cm 

Insulators  Used: 

EPEPEP 

EEEEEEEEEEEE 

EEEEEEEEEEEE 

i 

Xi  (cm) 

Ti  (K) 

Xi  (cm) 

Ti  (K) 

Xi  (cm) 

Ti  (K) 

1 

22.0 

8.38 

7.1875 

6.5875 

0.625 

4.25 

2 

23.8 

36.3 

11.406 

12.938 

8.125 

7.7375 

3 

24.8 

116.6 

15.625 

25.85 

7.9688 

12.369 

4 

29.4 

29.531 

71.094 

7.8125 

18.094 

5 

6.875 

100.31 

12.344 

29.912 

6 

5 

127.66 

26.094 

71.094 

7 

2.5 

143.13 

8.125 

105.94 

8 

2.5 

159.06 

5.3125 

135.47 

9 

4.6875 

188.59 

5 

165.94 

10 

5 

222.5 

5.625 

202.03 

11 

9.688 

12.9682 

We  can  see  immediately  that  the  addition  of  the  implicit  stress  constraint  results  in 
a  variable  cross-sectional  area  design  that  requires  over  50%  less  (normalized)  power  than 
the  design  found  by  Hilal  and  Eyssa  [23J.  This  savings  was  expected  because  the  newer 
formulation  allows  for  varying  the  number  of  heat  intercepts  and  the  mixing  of  insulator 
types.  A  similar  (65%)  savings  is  achieved  by  Kokkolaras  et  al.  |26j  in  optimizing  the 
constant  cross-sectional  area  formulation  of  Hilal  and  Boom  [22J,  but  with  the  additional 
categorical  design  variables. 

A  bit  surprising  is  the  slight  decrease  in  power  when  the  final  two  constraints  are  added 
(see  columns  3  and  4),  particularly  because  this  is  the  point  at  which  a  linearly  constrained 
problem  becomes  a  nonlinearly  constrained  one,  and  the  new  algorithm’s  filter  logic  is  ap¬ 
plied.  Recall  that  the  theory  ensures  convergence  to  a  first-order  stationary  point  for  the 
partial  model  (since  all  its  constraints  are  linear),  but  does  not  do  so  for  the  full  model.  In 
spite  of  this,  the  new  algorithm  still  finds  a  better  feasible  design.  It  is  indeed  possible  that 
the  FMGPS  algorithm  generates  a  different  sequence  and  simply  terminates  near  a  better 
local  minimizer. 
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Furthermore,  when  we  restart  the  process  with  this  solution  as  the  initial  point,  no 
movement  away  from  the  solution  occurs.  We  also  tried  starting  at  a  few  other  “common 
sense”  initial  designs,  including  both  (now  infeasible)  solutions  from  Table  2,  and  from 
a  design  with  10  evenly  spaced  intercepts  with  all  normal  epoxy  insulators;  however,  all  of 
these  produced  designs  that  required  more  power  than  the  full  model  design  given  in  Table  3, 
Thus  we  believe  our  solution  to  be  fairly  robust. 

Figure  5  illustrates  the  performance  of  the  FMGPS  algorithm  on  the  full  model,  where 
the  power  required  for  the  incumbent  best  design  is  plotted  versus  the  number  of  function 
evaluations.  The  lower  plot  is  a  magnification  of  the  upper  one.  The  “L” -shaped  plot  is  very 
typical  behavior  of  derivative-free  methods,  since  good  stopping  rules  for  these  methods 
are  difficult.  The  “stair  steps”  seen  in  the  right-hand  plot  indicate  varying  length  polling 
sequences. 

We  should  note  that  the  power  values  shown  on  the  vertical  axes  of  these  plots  do  not 
match  the  data  in  the  Table  3  because  they  represent  two  different  things.  The  objective 
function  is  to  minimize  power,  as  measured  in  Figure  5,  but  the  required  power  shown  in 
Table  3  is  normalized  (hence  the  (^)  notation),  so  as  to  allow  comparisons  with  the  results 
of  Hilal  and  Eyssa  |23j . 

Figure  6  depicts  the  progression  of  the  filter  during  the  run  of  the  full  model,  where  the 
plots  in  the  right  column  are  magnifications  of  those  on  the  left.  Each  of  the  three  rows 
represents  a  “snapshot”  taken  after  150,  200,  and  500  respective  function  evaluations  were 
performed.  Although  the  algorithm  terminated  after  more  than  9000  function  evaluations, 
changes  in  the  filter  after  500  function  evaluations  could  not  be  detected  within  the  resolution 
of  the  plot.  This  is  consistent  with  the  long  and  shallow  progression  of  the  best  objective 
function  value  seen  in  Figure  5.  Clearly,  better  stopping  rules  would  be  useful. 

In  the  filter  plots,  the  asterisks  represent  a  subset  of  the  best  feasible  points  found  up 
to  that  point,  while  the  “stair  step”  lines  represent  the  boundary  between  the  filtered  and 
unfiltered  points.  In  this  run,  the  nonlinear  constraints  were  scaled  by  dividing  each  by  its 
right-hand  side  and  then  subtracting  one  from  both  sides.  Thus  in  the  left  column  plots,  the 
choice  of  /rmax  =  1  represents  a  100%  constraint  violation. 

Table  3  also  shows  that  the  new  constraints  yield  significantly  different  insulator  configu¬ 
rations  than  that  of  Kokkolaras  et  al.  [26]  -  that  of  all  normal  cloth  fiberglass  epoxy.  This  is 
consistent  with  the  raw  data  we  used,  which  shows  epoxy  to  have  low  thermal  conductivity 
and  higher  resistance  to  stress  than  nylon  or  teflon. 

However,  some  of  the  other  materials  (including  nylon  and  teflon)  have  better  thermal 
contraction  properties  than  epoxy,  which  also  has  a  low  thermal  stress  threshold  and  which 
would  be  tested  by  any  thermal  contraction.  Modelling  thermal  stress  as  an  additional 
constraint  that  depends  on  thermal  contraction  would  be  an  interesting  extension  to  this 
problem  and  might  result  in  a  completely  different  material  configuration. 

The  results  shown  here  demonstrate  that,  although  the  FMGPS  algorithm  can  be  ex¬ 
pensive  when  applied  to  mixed  variable  problems,  it  successfully  generated  much-improved 
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Figure  5:  Iteration  History  for  the  Thermal  Insulation  System  Design  Problem 


designs  for  this  problem.  The  design  of  [23J  has  been  significantly  improved,  and  the  addition 
of  constraints  on  stress,  mass,  and  thermal  contraction  yields  a  more  realistic  feasible  design 
with  essentially  no  additional  power  required  over  that  of  [26J . 
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Figure  6:  Filter  Progression  for  the  Full  Model 
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